'Predator_Prey.bas
'program to plot the interaction between predator population (fox) and prey population (rabbit)

'set up the differential equations. fp=df/dt and rp=dr/dt
'fox populations grow with rabbits present but shrink with no rabbits
'rabbit populations decline when more foxes are around to eat them!
'the shape of the plot is very sensative to the exponent "w"

func fp(f,r)= .1*(f^w)*r-f
func rp(f,r)= 2*r-.05*f*r

'input initial conditions

input "exponent on fox population (.5 to 1) ";w
INPUT "initial number of foxes (10) "; f
INPUT "initial number of rabbits (100) "; r
input "duration of study (20) ";tf
n=400*tf
h = .0025
WINDOW -5, -20,n, 3*r

CLS

'draw graph axes

line 0,-20,0,3*r
line 0,0,n,0

'integrate the equtions

FOR q = 1 TO n
kfp1 = h * fp(f,r)
krp1 = h * rp(f,r)

kfp2 = h * fp(f + .5 * kfp1, r + .5 * krp1)
krp2 = h * rp(f + .5 * kfp1, r + .5 * krp1)

kfp3 = h * fp(f + .5 * kfp2, r + .5 * krp2)
krp3 = h * rp(f + .5 * kfp2, r + .5 * krp2)

kfp4 = h * fp(f + .5 * kfp3, r + .5 * krp3)
krp4 = h * rp(f + .5 * kfp3, r + .5 * krp3)

f = f + (kfp1 + 2 * kfp2 + 2 * kfp3 + kfp4) / 6
r = r + (krp1 + 2 * krp2 + 2 * krp3 + krp4) / 6

'plot the results

PSET q,f, 9: PSET q,r, 12
NEXT q
? "fox = blue"
? "rabbits = red"
Input "press enter to end ";bye
END